TIP24 Analysis
Zach Peagler
2024-09-23
Introduction
This analysis has been performed to accompany a tomato growth trial in the Hydroponic Lab at the KSU Field Station (coordinates?). The trial was conducted with the purpose of determining the effects of differing soil microbial inoculation timings on tomato plant photosynthesis and crop yield and quality. Photosynthesis metrics were taken with a Li-COR Li-600 and two PhotosynQ MultispeQ v2.0s, and include stomatal conductance (gsw), measured in (mol m-2 s-1), which is a measure of how well the leaf is able to conduct gas through its stomates, and the efficiency of photosystem II (PhiPS2), where photosystem II is a piece of machinery in the chloroplast essential to photosynthesis and PhiPS2 is the quantum yield in light calculated from fluorescence.
Sample size (n=48) with 4 groups, 12 replicates per group. The tomato plants were grown in a dutch bucket hydroponic system from April to October 2024. Fluorometric measurements were taken twice weekly, and fruit were harvested upong fully ripening and taken back to the lab on the KSU campus for analysis.
Fruit analysis included weighing the tomato for its mass in grams and visually assessing it for blossom end rot (BER), a nutritional disorder caused by calcium deficiency that renders the fruit unmarketable, as well as visual checks for fungus and cracking. Fungus commonly appears on parts of the tomato with blossom end rot or cracking. From there, fruit were analyzed with a penetrometer for ripeness, (1-kg) where kg is mapped from 0:1, and then with Fisher Brix Refractometers for sugar concentration in percent sugar content of the tomato’s juice.
Discarded fruit were then composted.
1. Setup
1.1 Load packages
DATA WRANGLING - Tidyverse
STATISTICS - lmerTest - Car - MASS - rstatix
GRAPHING - Showtext - Scico - ggpubr - ggridges
1.2 Graphical Setup
Load colors. Set the palette to one of the scientific color map palettes (scico) and create variants that contain the specific values for 2, 4, and 10 colors. This helps when working with factors, as it can be a pain to apply gradients otherwise. I also make the vector of colors =n+1 because the final value in scico palettes are often white, rendering them hard to see against a white plot background. (and we use white plot backgrounds because colored backgrounds are difficult to parse and expensive to print)
Load fonts. Google’s Open Sans for titles and Montserrat for body text.
Set shape selection from ggplot’s shape library.
a_palette <- "bilbao"
## it's a good idea to use n+1 for palettes as often the last color is white (invisible against the default background)
two_colors = scico(3, palette=a_palette)
four_colors = scico(5, palette=a_palette)
five_colors = scico(6, palette=a_palette)
true_two_col = scico(2, palette=a_palette)
ten_col = scico(10, palette=a_palette)
twelve_col = scico(13, palette=a_palette)
font_add_google("Open Sans", family = "open")
font_add_google("Montserrat", family = "mont")
showtext_auto()
four_shapes = c(15,16,17,23)2. Load Data
2.1 Li-600 Data
Load Li-600 data from file and convert strings to factors. Then change the treatments from numbers to the names of the groups. From here, filter the data for only those readings that have a leak percent of less than 10. QUESTION: What leak percent is acceptable to filter out? Rename the Date column to Date_ref to use it as a check that we formatted dates correctly. Format date to POSIXct in month-day-year format. Create a new plant factor that corresponds to each unique plant, then re-level the Treatment for consistency across analyses. Then we make a new column that is the log of gsw.
From here, outliers in the main target variables (stomatal conductance [gsw] and photosystem II efficiency [PhiPS2]) were identified. QUESTION: should I remove statistical outliers even if I assume the collection methods were not flawed.
Then, new dataframes were created out of the base Li-600 dataset to look at the mean and sd of gsw and PhiPS2, grouped by 1. treatment and date and 2. plant factor.
Another dataframe was created looking at mass and fruit count sums by treatment and plant.
As a final note, this data has something wrong with it, as there are ~100 gsw values <0, which should be impossible. We’ve contacted Li-COR for clarification and hopefully there’s just a correction we can apply to the data. If not, then well uhhhh…-o.o-
## Load and clean Li-600 data
Li_data_file <- "C:/Github/Portfolio/R Shiny Apps/01 Tomato Inoculants 2024/TIP24_LI600.csv"
Li_data <- read.csv(Li_data_file, stringsAsFactors = T) %>%
mutate(Treatment = case_when(
Row==1~"Control",
Row==2~"Transplantation",
Row==3~"Germination",
Row==4~"Germ+Trans",
TRUE~NA)) %>%
filter(leak_pct<10, gsw>0) %>%
rename(Date_ref = Date) %>%
mutate(Date = parse_date_time(Date_ref, orders = "mdy"),
plant_fac = as.factor(paste(Row, Pot)),
Treatment = factor(Treatment, levels=c("Control", "Germination",
"Transplantation", "Germ+Trans")),
loggsw = log(gsw)
)
# Identify outliers of target variables
## Get gsw and PhiPS2 quantiles and interquartile range (IQR)
Qgsw <- quantile(Li_data$gsw, probs=c(.25, .75), na.rm=FALSE)
Qps2 <- quantile(Li_data$PhiPS2, probs=c(.25, .75), na.rm=FALSE)
## Get interquartile range
iqr_gsw <- IQR(Li_data$gsw)
iqr_ps2 <- IQR(Li_data$PhiPS2)
## Make new dataframes without outliers
Li_data_no_outliers <- subset(Li_data, Li_data$gsw > (Qgsw[1]-1.5*iqr_gsw) &
Li_data$gsw < (Qgsw[2]+1.5*iqr_gsw) &
Li_data$PhiPS2 > (Qps2[1]-1.5*iqr_ps2) &
Li_data$PhiPS2 < (Qps2[2]+1.5*iqr_ps2)
)
## Make new dataframes without extreme outliers
Li_data_no_extremes <- subset(Li_data, Li_data$gsw > (Qgsw[1]-3.0*iqr_gsw) &
Li_data$gsw < (Qgsw[2]+3.0*iqr_gsw) &
Li_data$PhiPS2 > (Qps2[1]-3.0*iqr_ps2) &
Li_data$PhiPS2 < (Qps2[2]+3.0*iqr_ps2)
)
# New dataframes of Li_data mean and sd for stomatal conductance (gsw) and Photosystem II efficiency (PhiPS2), grouped by treatment and plant
Li_data_stats <- Li_data %>%
group_by(Treatment) %>%
summarise_at(vars(gsw, PhiPS2), list(mean=mean, sd=sd))
Li_data_stats_byplant <- Li_data %>%
group_by(Treatment, plant_fac) %>%
summarise_at(vars(gsw, PhiPS2), list(mean=mean, sd=sd))
## Print mean and sd for Li_data by treatment
print(Li_data_stats)## # A tibble: 4 × 5
## Treatment gsw_mean PhiPS2_mean gsw_sd PhiPS2_sd
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Control 0.696 0.703 0.550 0.0415
## 2 Germination 0.523 0.716 0.511 0.0525
## 3 Transplantation 0.562 0.724 0.486 0.0637
## 4 Germ+Trans 0.404 0.730 0.336 0.0568#print(Li_data_stats_byplant)2.2 Fruit Data
Load fruit data from set file name and location, converting strings to factors and renaming Treatment’s group names. Then setup “fruit” as a way to sum the total fruit count, create references for our dates then parse them to POSIXct in month-day-year format, then grab just the days and calculate the difference in days from harvest to analysis. (NOTE: this is susceptible to month rollover errors, so be sure not to make any of those mistakes) Make a plant factor as a unique ID for each plant, relevel the treatment for consistency across analyses, create mass bins (n=10) to look at BER occurrence across mass groups, and convert BER, fungus, and cracking to factors. Then, we map the penetrometer data from 0:1 and subtract it from 1.0, effectively giving us the ripeness. (this is because higher penetrometer values correspond to lower ripeness)
Then create a copy of the data that doesn’t have BER and another copy of just those with BER.
Make a pivot table of the means and sd of the tomato mass (g) and sugar (%), grouped by treatment.
Make another pivot table to get the sum of fruit count, mass, and BER, fungus, and cracking occurrence by treatment group. Then calculate the probability of BER, fungus, and cracking for each group. Do that again grouping the fruit into mass bins.
## Load and clean TIP24 fruit data
Fl_data_file <- "C:/Github/Portfolio/R Shiny Apps/01 Tomato Inoculants 2024/TIP24_Fruit.csv"
Fl_data <- read.csv(Fl_data_file, stringsAsFactors = T) %>%
mutate(Treatment = case_when(
row==1~"Control",
row==2~"Transplantation",
row==3~"Germination",
row==4~"Germ+Trans",
TRUE~NA),
fruit = 1,
date_analysis_ref = date_analysis,
date_harvest_ref = date_harvest,
date_analysis = parse_date_time(date_analysis_ref, orders = "mdy"),
date_harvest = parse_date_time(date_harvest_ref, orders="mdy"),
d_harvest = format(date_analysis, format="%d"),
d_analysis = format(date_harvest, format="%d"),
d_diff = abs(as.integer(d_analysis)-as.integer(d_harvest)),
plant_fac = as.factor(paste(row, plant)),
Treatment = factor(Treatment, levels=c("Control","Germination",
"Transplantation",
"Germ+Trans")),
mass_bin = cut(mass, breaks=10),
BER_fac = as.factor(BER),
fungus_fac = as.factor(fungus),
cracking_fac = as.factor(cracking)
)
## Filter to only rows where there is no BER and sugar is greater than 0
Fl_data_no_BER <- Fl_data %>%
filter(BER==0 & sugar_avg > 0) %>%
mutate(pen_mapped = penetrometer/max(penetrometer),
ripeness = abs(1-pen_mapped)
)
# Identify outliers in target variables
## Get mass and sugar quantiles and interquartile range (IQR)
Qmass <- quantile(Fl_data_no_BER$mass, probs=c(.25, .75), na.rm=FALSE)
Qsug <- quantile(Fl_data_no_BER$sugar_avg, probs=c(.25, .75), na.rm=FALSE)
## Get interquartile range
iqr_mass <- IQR(Fl_data_no_BER$mass)
iqr_sug <- IQR(Fl_data_no_BER$sugar_avg)
## Make new dataframes without outliers
Fl_data_no_outliers <- subset(Fl_data_no_BER, mass > (Qmass[1]-1.5*iqr_mass) &
mass < (Qmass[2]+1.5*iqr_mass) &
sugar_avg > (Qsug[1]-1.5*iqr_sug) &
sugar_avg < (Qsug[2]+1.5*iqr_sug)
)
## Make new dataframes without extremes
Fl_data_no_extremes <- subset(Fl_data_no_BER, mass > (Qmass[1]-3.0*iqr_mass) &
mass < (Qmass[2]+3.0*iqr_mass) &
sugar_avg > (Qsug[1]-3.0*iqr_sug) &
sugar_avg < (Qsug[2]+3.0*iqr_sug)
)
# New dataframes
## New dataframes, filtering to only those where BER is equal to 1
Fl_data_BER <- Fl_data %>%
filter(BER==1)
## Make summary dataframe with mean and sd of mass and sugar by treatment.
Fl_data_means <- Fl_data_no_BER %>%
group_by(Treatment) %>%
summarise_at(vars(mass, sugar_avg), list(mean=mean, sd=sd))
## Make summary dataframe with mean and sd of mass and sugar by plant.
Fl_means_byplant <- Fl_data_no_BER %>%
group_by(Treatment, plant_fac) %>%
summarise_at(vars(mass, sugar_avg), list(mean=mean, sd=sd))
## Sum the fruit, BER, fungus, cracking, and mass by treatment group, then get probs of each factor.
Fl_data_summary <- Fl_data %>%
group_by(Treatment) %>%
summarise_at(vars(fruit, mass, BER, fungus, cracking, mass),
list(sum=sum)) %>%
mutate(pBER = round(BER_sum/fruit_sum, 4),
pfungus = round(fungus_sum/fruit_sum, 4),
pcracking = round(cracking_sum/fruit_sum, 4)
)
Fl_data_summary_date <- Fl_data %>%
group_by(Treatment, date_analysis) %>%
summarise_at(vars(fruit, mass, BER, fungus, cracking, mass),
list(sum=sum)) %>%
mutate(pBER = round(BER_sum/fruit_sum, 4),
pfungus = round(fungus_sum/fruit_sum, 4),
pcracking = round(cracking_sum/fruit_sum, 4)
)
## Sum the fruit, BER, fungus, and cracking by mass bin, then get probs of each factor.
Fl_data_mb <- Fl_data %>%
group_by(Treatment, mass_bin) %>%
summarise_at(vars(fruit, mass, BER, fungus, cracking), list(sum=sum)) %>%
mutate(pBER = round(BER_sum/fruit_sum, 4)*100,
pfungus = round(fungus_sum/fruit_sum, 4)*100,
pcracking = round(cracking_sum/fruit_sum, 4)*100
)
## Sum the fruit, BER, fungus, and cracking by plant, then get probs of each factor.
Fl_sum_byplant <- Fl_data %>%
group_by(Treatment, plant) %>%
summarise_at(vars(fruit, mass, BER, fungus, cracking), list(sum=sum)) %>%
mutate(pBER = round(BER_sum/fruit_sum, 4)*100,
pfungus = round(fungus_sum/fruit_sum, 4)*100,
pcracking = round(cracking_sum/fruit_sum, 4)*100
)
## Print fruit data means
print(Fl_data_means)## # A tibble: 4 × 5
## Treatment mass_mean sugar_avg_mean mass_sd sugar_avg_sd
## <fct> <dbl> <dbl> <dbl> <dbl>
## 1 Control 120. 6.43 68.6 1.54
## 2 Germination 130. 6.71 81.4 1.47
## 3 Transplantation 131. 6.66 76.3 1.26
## 4 Germ+Trans 106. 6.98 75.1 1.71## Print fruit data means by plant
#print(Fl_means_byplant)
## Print fruit summary
print(Fl_data_summary)## # A tibble: 4 × 9
## Treatment fruit_sum mass_sum BER_sum fungus_sum cracking_sum pBER pfungus
## <fct> <dbl> <dbl> <int> <int> <int> <dbl> <dbl>
## 1 Control 197 22352. 25 56 62 0.127 0.284
## 2 Germination 263 21677. 34 38 50 0.129 0.144
## 3 Transplanta… 361 29968. 29 51 98 0.0803 0.141
## 4 Germ+Trans 204 19191. 22 32 53 0.108 0.157
## # ℹ 1 more variable: pcracking <dbl>## Print fruit summary by mass bin
print(Fl_data_mb)## # A tibble: 37 × 10
## # Groups: Treatment [4]
## Treatment mass_bin fruit_sum mass_sum BER_sum fungus_sum cracking_sum pBER
## <fct> <fct> <dbl> <dbl> <int> <int> <int> <dbl>
## 1 Control (3.59,4… 28 904. 5 8 5 17.9
## 2 Control (45.2,8… 51 3412. 6 12 17 11.8
## 3 Control (86.5,1… 57 5967. 11 20 20 19.3
## 4 Control (128,16… 24 3573. 1 5 7 4.17
## 5 Control (169,21… 16 3050. 1 8 8 6.25
## 6 Control (210,25… 12 2807. 0 3 3 0
## 7 Control (252,29… 5 1344. 0 0 2 0
## 8 Control (293,33… 3 939. 1 0 0 33.3
## 9 Control (334,37… 1 356. 0 0 0 0
## 10 Germination (3.59,4… 118 3265. 24 17 14 20.3
## # ℹ 27 more rows
## # ℹ 2 more variables: pfungus <dbl>, pcracking <dbl>## Print fruit summary by plant
#print(Fl_sum_byplant)3. Analyses
3.1 Li-600 Analysis
3.1.1 Exploratory Li-600 Analysis
Check for interactions and get a feel for the data. Use Q-Q plots and Shapiro-Wilk tests to check gsw and PhiPS2 for normality.
QUESTION: How do I change the xlim when x is a non-numeric factor.
# check for normality
shapiro.test(Li_data$PhiPS2)##
## Shapiro-Wilk normality test
##
## data: Li_data$PhiPS2
## W = 0.75205, p-value < 2.2e-16shapiro.test(Li_data$gsw)##
## Shapiro-Wilk normality test
##
## data: Li_data$gsw
## W = 0.7817, p-value < 2.2e-16## gsw distributions
gsw_n <- fitdistr(Li_data$gsw, "normal")
gsw_ln <- fitdistr(Li_data$gsw, "lognormal")
gsw_g <- fitdistr(Li_data$gsw, "gamma")
gsw_exp <- fitdistr(Li_data$gsw, "exponential")
gsw_seq <- seq(min(Li_data$gsw), max(Li_data$gsw), by=0.001)
## PDFs
gsw_pdf_n <- dnorm(gsw_seq, mean=gsw_n$estimate[1],
sd=gsw_n$estimate[2])
gsw_pdf_ln <- dlnorm(gsw_seq, meanlog=gsw_ln$estimate[1],
sdlog=gsw_ln$estimate[2])
gsw_pdf_g <- dgamma(gsw_seq, shape=gsw_g$estimate[1],
rate=gsw_g$estimate[2])
gsw_pdf_exp <- dexp(gsw_seq, rate=gsw_exp$estimate)
gsw_pdf <- density(Li_data$gsw, n=length(gsw_pdf_ln))
gsw_pdf_df <- as.data.frame(gsw_seq) %>%
mutate(pdf_ln = gsw_pdf_ln,
pdf_n = gsw_pdf_n,
pdf_g = gsw_pdf_g,
pdf_exp = gsw_pdf_exp)
## CDFs
gsw_cdf_n <- pnorm(gsw_seq, mean=gsw_n$estimate[1],
sd=gsw_n$estimate[2])
gsw_cdf_ln <- plnorm(gsw_seq, meanlog=gsw_ln$estimate[1],
sdlog=gsw_ln$estimate[2])
gsw_cdf_g <- pgamma(gsw_seq, shape=gsw_g$estimate[1],
rate=gsw_g$estimate[2])
gsw_cdf_exp <- pexp(gsw_seq, rate=gsw_exp$estimate)
gsw_cdf <- ecdf(Li_data$gsw)(gsw_seq)
gsw_cdf_df <- as.data.frame(gsw_seq) %>%
mutate(cdf_ln = gsw_cdf_ln,
cdf_n = gsw_cdf_n,
cdf_g = gsw_cdf_g,
cdf_exp = gsw_cdf_exp,
cdf = gsw_cdf)
ggplot(gsw_pdf_df)+
geom_point(aes(x=gsw_seq, y=pdf_n, color="Normal"))+
geom_point(aes(x=gsw_seq, y=pdf_ln, color="Lognormal"))+
geom_point(aes(x=gsw_seq, y=pdf_g, color="Gamma"))+
geom_point(aes(x=gsw_seq, y=pdf_exp, color="Exponential"))+
geom_point(aes(x=gsw_pdf$x, y=gsw_pdf$y, color="Data"))+
theme_bw()+
scale_color_manual(values=five_colors)+
labs(title="PDF for gsw data and possible distributions")+
ylab("PDF")+
xlab("gsw")+
theme(
text = element_text(size=24, family="mont"),
legend.position="inside",
legend.position.inside = c(.85,.6),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
ggplot(gsw_cdf_df)+
geom_point(aes(x=gsw_seq, y=cdf_n, color="Normal"))+
geom_point(aes(x=gsw_seq, y=cdf_ln, color="Lognormal"))+
geom_point(aes(x=gsw_seq, y=cdf_g, color="Gamma"))+
geom_point(aes(x=gsw_seq, y=cdf_exp, color="Exponential"))+
geom_point(aes(x=gsw_seq, y=cdf, color="Data"))+
theme_bw()+
scale_color_manual(values=five_colors)+
labs(title="CDF for gsw data and possible distributions")+
ylab("CDF")+
xlab("gsw")+
theme(
text = element_text(size=24, family="mont"),
legend.position="inside",
legend.position.inside = c(.85,.6),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
### test gsw for gamma and plnorm distributions using ks test
gsw_gks<- ks.test(Li_data$gsw, "pgamma",shape=gsw_g$estimate[1],
rate=gsw_g$estimate[2])
ks.test(Li_data$gsw, "plnorm",meanlog=gsw_ln$estimate[1],
sdlog=gsw_ln$estimate[2])##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: Li_data$gsw
## D = 0.088853, p-value = 0.0002042
## alternative hypothesis: two-sided## PhiPS2 distribution
PhiPS2_n <- fitdistr(Li_data$PhiPS2, "normal")
PhiPS2_ln <- fitdistr(Li_data$PhiPS2, "lognormal")
PhiPS2_g <- fitdistr(Li_data$PhiPS2, "gamma")
PhiPS2_exp <- fitdistr(Li_data$PhiPS2, "exponential")
PhiPS2_seq <- seq(min(Li_data$PhiPS2), max(Li_data$PhiPS2), by=0.001)
## PDFs
PhiPS2_pdf_n <- dnorm(PhiPS2_seq, mean=PhiPS2_n$estimate[1],
sd=PhiPS2_n$estimate[2])
PhiPS2_pdf_ln <- dlnorm(PhiPS2_seq, meanlog=PhiPS2_ln$estimate[1],
sdlog=PhiPS2_ln$estimate[2])
PhiPS2_pdf_g <- dgamma(PhiPS2_seq, shape=PhiPS2_g$estimate[1],
rate=PhiPS2_g$estimate[2])
PhiPS2_pdf_exp <- dexp(PhiPS2_seq, rate=PhiPS2_exp$estimate)
PhiPS2_pdf <- density(Li_data$PhiPS2, n=length(PhiPS2_pdf_ln))
PhiPS2_pdf_df <- as.data.frame(PhiPS2_seq) %>%
mutate(pdf_ln = PhiPS2_pdf_ln,
pdf_n = PhiPS2_pdf_n,
pdf_g = PhiPS2_pdf_g,
pdf_exp = PhiPS2_pdf_exp)
## CDFs
PhiPS2_cdf_n <- pnorm(PhiPS2_seq, mean=PhiPS2_n$estimate[1],
sd=PhiPS2_n$estimate[2])
PhiPS2_cdf_ln <- plnorm(PhiPS2_seq, meanlog=PhiPS2_ln$estimate[1],
sdlog=PhiPS2_ln$estimate[2])
PhiPS2_cdf_g <- pgamma(PhiPS2_seq, shape=PhiPS2_g$estimate[1],
rate=PhiPS2_g$estimate[2])
PhiPS2_cdf_exp <- pexp(PhiPS2_seq, rate=PhiPS2_exp$estimate)
PhiPS2_cdf <- ecdf(Li_data$PhiPS2)(PhiPS2_seq)
PhiPS2_cdf_df <- as.data.frame(PhiPS2_seq) %>%
mutate(cdf_ln = PhiPS2_cdf_ln,
cdf_n = PhiPS2_cdf_n,
cdf_g = PhiPS2_cdf_g,
cdf_exp = PhiPS2_cdf_exp,
cdf = PhiPS2_cdf)
## PhiPS2 PDF plot
ggplot(PhiPS2_pdf_df)+
geom_point(aes(x=PhiPS2_seq, y=pdf_n, color="Normal"))+
geom_point(aes(x=PhiPS2_seq, y=pdf_ln, color="Lognormal"))+
geom_point(aes(x=PhiPS2_seq, y=pdf_g, color="Gamma"))+
#geom_point(aes(x=PhiPS2_seq, y=pdf_exp, color="Exponential"))+
geom_point(aes(x=PhiPS2_pdf$x, y=PhiPS2_pdf$y, color="Data"))+
theme_bw()+
scale_color_manual(values=five_colors)+
guides(color=guide_legend(title="Distribution"))+
labs(title="PDF for PhiPS2 data and possible distributions")+
ylab("PDF")+
xlab("PhiPS2")+
theme(
text = element_text(size=24, family="mont"),
legend.position="inside",
legend.position.inside = c(.85,.6),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
## PhiPS2 CDF plot
ggplot(PhiPS2_cdf_df)+
geom_point(aes(x=PhiPS2_seq, y=cdf_n, color="Normal"))+
geom_point(aes(x=PhiPS2_seq, y=cdf_ln, color="Lognormal"))+
geom_point(aes(x=PhiPS2_seq, y=cdf_g, color="Gamma"))+
geom_point(aes(x=PhiPS2_seq, y=cdf, color="Data"))+
theme_bw()+
scale_color_manual(values=five_colors)+
guides(color=guide_legend(title="Distribution"))+
labs(title="CDF for PhiPS2 data and possible distributions")+
ylab("CDF")+
xlab("PhiPS2")+
theme(
text = element_text(size=24, family="mont"),
legend.position="inside",
legend.position.inside = c(.85,.6),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
### test PhiPS2 for gamma , norm, and lnorm distributions using ks test
ks.test(Li_data$PhiPS2, "pgamma",shape=PhiPS2_g$estimate[1],
rate=PhiPS2_g$estimate[2])## Warning in ks.test.default(Li_data$PhiPS2, "pgamma", shape =
## PhiPS2_g$estimate[1], : ties should not be present for the Kolmogorov-Smirnov
## test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: Li_data$PhiPS2
## D = 0.13486, p-value = 1.279e-09
## alternative hypothesis: two-sidedks.test(Li_data$PhiPS2, "plnorm",meanlog=PhiPS2_ln$estimate[1],
sdlog=PhiPS2_ln$estimate[2])## Warning in ks.test.default(Li_data$PhiPS2, "plnorm", meanlog =
## PhiPS2_ln$estimate[1], : ties should not be present for the Kolmogorov-Smirnov
## test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: Li_data$PhiPS2
## D = 0.11891, p-value = 1.424e-07
## alternative hypothesis: two-sidedks.test(Li_data$PhiPS2, "pnorm",mean=PhiPS2_n$estimate[1],
sd=PhiPS2_n$estimate[2])## Warning in ks.test.default(Li_data$PhiPS2, "pnorm", mean =
## PhiPS2_n$estimate[1], : ties should not be present for the Kolmogorov-Smirnov
## test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: Li_data$PhiPS2
## D = 0.14124, p-value = 1.646e-10
## alternative hypothesis: two-sided# check for heterodescasticity
leveneTest(Li_data$gsw~Li_data$Treatment)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.8428 0.03717 *
## 578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(Li_data$PhiPS2~Li_data$Treatment)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 3.6563 0.01242 *
## 578
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Stomatal Conductance density by treatment
ggplot(data=Li_data, aes(x=gsw, y = Treatment, color=Treatment, fill=Treatment))+
geom_density_ridges()+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
#xlim(-0.25, 1.5)+
ylab("Density")+
xlab("Stomatal Conductance (mol m-2 s-1)")+
theme_minimal()+
labs(
title=str_wrap("Stomatal Conductance Density for Inoculation Treatments in Tomato", 40)
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight = 0.5),
axis.title.y = element_text(size=0),
axis.text = element_text(size=24, family="mont"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Picking joint bandwidth of 0.0964
## gsw Density Graph by plant
ggplot(data=Li_data, aes(x=gsw, y = plant_fac, color=Treatment, fill=Treatment))+
geom_density_ridges()+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
#xlim(-.025, 1.5)+
ylab("Density")+
xlab("Stomatal Conductance (mol m-2 s-1)")+
theme_minimal()+
labs(
title=str_wrap("Stomatal Conductance Density for Different Tomato Plants by Treatment", 40)
)+
theme(
legend.position="right",
text = element_text(size=24, family="mont", lineheight = 0.5),
axis.title.y = element_text(size=0),
axis.text.x = element_text(size=24, family="mont"),
axis.text = element_text(size=12, family="mont"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Picking joint bandwidth of 0.141
## PhiPS2 Density Graph by treatment
ggplot(data=Li_data, aes(x=PhiPS2, y = Treatment, color=Treatment, fill=Treatment))+
geom_density_ridges()+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
#xlim(.6,.8)+
ylab("Density")+
xlab("Efficiency of Photosystem II")+
theme_minimal()+
labs(
title=str_wrap("Photosystem II Efficiency Density for Inoculation Treatments in Tomato", 40)
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight = 0.5),
axis.title.y = element_text(size=0),
axis.text = element_text(size=24, family="mont"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Picking joint bandwidth of 0.0112
## PhiPS2 Density Graph by plant
ggplot(data=Li_data, aes(x=PhiPS2, y = plant_fac, color=Treatment, fill=Treatment))+
geom_density_ridges()+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
xlim(.6,.8)+
ylab("Density")+
xlab("Stomatal Conductance (mol m-2 s-1)")+
theme_minimal()+
labs(
title=str_wrap("Photosystem II Efficiency Density for Different Tomato Plants by Treatment", 40)
)+
theme(
legend.position="right",
text = element_text(size=24, family="mont", lineheight = 0.5),
axis.title.y = element_text(size=0),
axis.text.x = element_text(size=24, family="mont"),
axis.text = element_text(size=12, family="mont"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Picking joint bandwidth of 0.0148## Warning: Removed 24 rows containing non-finite outside the scale range
## (`stat_density_ridges()`).
# check for correlation between PhiPS2 and Qamb
ggplot(data = Li_data, aes(x=Qamb, y = PhiPS2,
shape=Treatment, fill=Treatment, color = Treatment)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
scale_fill_manual(values=four_colors)+
ylim(.6,.8)+
theme_minimal()+
ylab("PhiPS2")+
xlab("Ambient Light")+
labs(
title=str_wrap("Photosystem II Efficiency By Ambient Light Across Inoculation Treatment in Tomato", 50)
)+
theme(
text = element_text(size=24, family="mont"),
legend.title = element_text(size=24, family="mont", face="bold"),
legend.text = element_text(size=20, family="mont"),
legend.position="inside",
legend.title.position = "top",
legend.position.inside=c(0.15,0.2),
legend.key.height = unit(.3, "cm"),
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Plot PhiPS2 against Qamb, facet wrapped by treatment and colored by Temperature. Because Qamb and Temperature are almost certainly correlated, this may not be very useful.
ggplot(data = Li_data, aes(x=Qamb, y = PhiPS2, color = Tleaf)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_scico(begin=0.9, end=0, palette=a_palette)+
facet_wrap(~Treatment)+
ylim(.6,.8)+
theme_minimal()+
ylab("PhiPS2")+
xlab("Ambient Light")+
labs(
title=str_wrap("Photosystem II Efficiency By Ambient Light Across Temperature By Inoculation Treatment in Tomato", 50)
)+
guides(color=guide_colorbar(title=str_wrap("Leaf Temp (C)", 8)))+
theme(
text = element_text(size=24, family="mont"),
legend.position="right",
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
3.1.2 Stomatal Conductance Analysis
Do a generalized linear mixed model with stomatal conductance (gsw) as response, Treatment as predictor, relative humidity (rh_s) as fixed effect, and plant as random effect, with an inverse link function because gsw has a gamma distribution.
#GLMM with Treatment as explanatory
#summary(gsw_glmm <- (glmer(gsw ~ Treatment + rh_s + (1|plant_fac),
# data=Li_data, family=Gamma(link="log"))))
# GSW by date across relative humidity
ggplot(data = Li_data, aes(x=as.factor(Date), y = gsw,
color = rh_s)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_scico(begin=0.9, end=0, palette=a_palette)+
scale_x_discrete(guide=guide_axis(check.overlap=TRUE))+
ylim(0, 1.5)+
theme_minimal()+
ylab("Stomatal Conductance (mol m-2 s-1)")+
xlab("Date")+
labs(
title=str_wrap("Stomatal Conductance By Date Across Relative Humidity in Tomato", 50)
)+
guides(color=guide_colorbar(title=str_wrap("Relative Humidity", 10)))+
theme(
text = element_text(size=24, family="mont"),
legend.title = element_text(size=24, family="mont", face="bold"),
legend.text = element_text(size=20, family="mont"),
legend.position="inside",
legend.title.position = "top",
legend.position.inside=c(0.8,0.8),
legend.key.height = unit(.3, "cm"),
legend.background = element_rect(color=two_colors[2], fill=NA,
linewidth=.5, linetype = 2),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Warning: Removed 31 rows containing missing values or values outside the scale range
## (`geom_point()`).
# gsw by relative humidity across treatment groups
ggplot(data = Li_data, aes(x=rh_s, y = log(gsw),
shape=Treatment, fill=Treatment, color = Treatment)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
scale_fill_manual(values=four_colors)+
#ylim(0,1)+
theme_minimal()+
ylab("Log of Stomatal Conductance (mol m-2 s-1)")+
xlab("Relative Humidity")+
labs(
title=str_wrap("Log Adjusted Stomatal Conductance By Relative Humidity Across Inoculation Treatment in Tomato", 50)
)+
theme(
text = element_text(size=24, family="mont"),
legend.title = element_text(size=24, family="mont", face="bold"),
legend.text = element_text(size=20, family="mont"),
legend.position="inside",
legend.title.position = "top",
legend.position.inside=c(0.85,0.15),
legend.key.height = unit(.3, "cm"),
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
# The following boxplots are probably not the best way to evaluate gsw.
## gsw by Treatment boxplot
ggplot(data = Li_data, aes(x= Treatment, y = gsw, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.25)+
geom_violin(alpha=0.5, width=1)+
geom_jitter( width=.2, height=0)+
stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
stat_compare_means(label.x=3, label.y=3.5, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("Stomatal Conductance Across Inoculation Treatments in Tomato", 40)
) +
ylab("Stomatal Conductance (mol m-2 s-1)")+
theme_minimal()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## gsw by Treatment boxplot, NO OUTLIERS
ggplot(data = Li_data_no_outliers, aes(x= Treatment, y = gsw, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.25)+
geom_violin(alpha=0.5, width=1)+
geom_jitter( width=.2, height=0)+
stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
stat_compare_means(label.x=3, label.y=1.5, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("Stomatal Conductance Across Inoculation Treatments in Tomato No Outliers", 40)
) +
ylab("Stomatal Conductance (mol m-2 s-1)")+
theme_minimal()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
# Mean GSW per plant by Treatment boxplot
ggplot(data = Li_data_stats_byplant, aes(x= Treatment, y = gsw_mean, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.1)+
geom_violin(alpha=0.5, width=1.5)+
geom_jitter( width=.2, height=0)+
stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans"), c("Germination", "Transplantation"), c("Germination", "Germ+Trans"), c("Transplantation", "Germ+Trans")), size=8, family="mont")+
stat_compare_means(label.x=3, label.y=.8, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("Mean Stomatal Conductance Per Plant Across Inoculation Treatments in Tomato", 50)
) +
ylab("Stomatal Conductance (mol m-2 s-1)")+
theme_minimal()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Warning: `position_dodge()` requires non-overlapping x intervals.
3.1.3 PhiPS2 Analysis
# linear mixed model with PhiPS2 as response, Treatment as predictor (explanatory), light (Qamb) as fixed effect, and plant as random effect
#summary(model_PhiPS2 <- (lmer(
# PhiPS2 ~ Treatment + Qamb + (1|plant_fac),
# data = Li_data)))
#summary(model_PhiPS2_noT <- (lmer(
# PhiPS2 ~ Qamb + (1|plant_fac),
# data = Li_data)))
## PhiPS2 by date, grouped by treatment
ggplot(data = Li_data, aes(x=as.factor(Date), y = PhiPS2, color = Treatment,
fill=Treatment, shape=Treatment)) +
geom_jitter(width=1, size=2, alpha=1)+
scale_color_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
scale_fill_manual(values=four_colors)+
scale_x_discrete(guide=guide_axis(check.overlap=TRUE))+
ylim(.6,.8)+
facet_wrap(~Treatment)+
theme_minimal()+
ylab("PhiPS2")+
xlab("Date")+
labs(
title=str_wrap("PhiPS2 By Date Across Microbial Treatments in Tomato", 60)
)+
theme(
text = element_text(size=24, family="mont"),
legend.position="none",
axis.title = element_text(size=24, family = "mont", face= "bold"),
axis.text.x = element_text(size=12, family = "mont", angle=10),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
# PhiPS2 by leaf across temperature
ggplot(data = Li_data, aes(x=Date, y = PhiPS2, color = Tleaf)) +
geom_jitter(width=1, size=2, alpha=1, aes(x=as.factor(Date)))+
scale_color_scico(begin=0.9, end=0, palette=a_palette)+
scale_x_discrete(guide=guide_axis(check.overlap=TRUE))+
ylim(.6,.8)+
theme_minimal()+
guides(color=guide_legend(title="Leaf Temp (C)"))+
ylab("PhiPS2")+
xlab("Date")+
labs(
title=str_wrap("PhiPS2 By Date Across Temperature in Tomato", 50)
)+
theme(
text = element_text(size=24, family="mont"),
legend.title = element_text(size=24, family="mont", face="bold"),
legend.text = element_text(size=20, family="mont"),
legend.position="inside",
legend.title.position = "top",
legend.position.inside=c(0.75,0.2),
legend.key.height = unit(.3, "cm"),
legend.background = element_rect(color=two_colors[2], fill=NA, linewidth=.5, linetype = 2),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)## Warning: Removed 24 rows containing missing values or values outside the scale range
## (`geom_point()`).
# The following boxplots are probably not the best way to evaluate gsw.
## PhiPS2 by Treatment boxplot
ggplot(data = Li_data, aes(x= Treatment, y = PhiPS2, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.25)+
geom_violin(alpha=0.5, width=1)+
geom_jitter( width=.2, height=0)+
stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
stat_compare_means(label.x=3.5, label.y=1.05, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("Efficiency of Photosystem II Across Inoculation Treatments in Tomato", 40)
) +
ylab("PhiPS2")+
theme_minimal()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## PhiPS2 by Treatment boxplot NO OUTLIERS
ggplot(data = Li_data_no_outliers, aes(x= Treatment, y = PhiPS2, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.25)+
geom_violin(alpha=0.5, width=1)+
geom_jitter( width=.2, height=0)+
stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
stat_compare_means(label.x=3.25, label.y=.85, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("Efficiency of Photosystem II Across Inoculation Treatments in Tomato No Outliers", 40)
) +
ylab("PhiPS2")+
theme_minimal()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## Mean PhiPS2 per plant by Treatment boxplot
ggplot(data = Li_data_stats_byplant, aes(x= Treatment, y = PhiPS2_mean, fill=Treatment, color=Treatment)) +
geom_boxplot(alpha=.5, width=0.1)+
geom_violin(alpha=0.5, width=1)+
geom_jitter( width=.2, height=0)+
stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans"), c("Germination", "Transplantation"), c("Germination", "Germ+Trans"), c("Transplantation", "Germ+Trans")), size=8, family="mont")+
stat_compare_means(label.x=.5, label.y=.8, size=8,family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
labs(
title=str_wrap("Mean Photosystem II Efficiency Per Plant Across Inoculation Treatments in Tomato", 50)
) +
ylab("PhiPS2")+
theme_minimal()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
3.2 Fruit Analysis
3.2.1 Exploratory Fruit Analysis
# check for normality
## use shapiro wilk tests
shapiro.test(Fl_data_no_BER$mass)##
## Shapiro-Wilk normality test
##
## data: Fl_data_no_BER$mass
## W = 0.91001, p-value < 2.2e-16shapiro.test(Fl_data_no_BER$sugar_avg)##
## Shapiro-Wilk normality test
##
## data: Fl_data_no_BER$sugar_avg
## W = 0.98414, p-value = 2.486e-06## mass distribution
mass_n <- fitdistr(Fl_data$mass, "normal")
mass_ln <- fitdistr(Fl_data$mass, "lognormal")
mass_g <- fitdistr(Fl_data$mass, "gamma")
mass_seq <- seq(min(Fl_data$mass), max(Fl_data$mass), by=0.1)
mass_pdf_n <- dnorm(mass_seq, mean=mass_n$estimate[1],
sd=mass_n$estimate[2])
mass_pdf_ln <- dlnorm(mass_seq, meanlog=mass_ln$estimate[1],
sdlog=mass_ln$estimate[2])
mass_pdf_g <- dgamma(mass_seq, shape=mass_g$estimate[1],
rate=mass_g$estimate[2])
mass_pdf <- density(Fl_data$mass, n=length(mass_pdf_ln))
mass_pdf_df <- as.data.frame(mass_seq) %>%
mutate(pdf_ln = mass_pdf_ln,
pdf_n = mass_pdf_n,
pdf_g = mass_pdf_g)
## CDFs
mass_cdf_n <- pnorm(mass_seq, mean=mass_n$estimate[1],
sd=mass_n$estimate[2])
mass_cdf_ln <- plnorm(mass_seq, meanlog=mass_ln$estimate[1],
sdlog=mass_ln$estimate[2])
mass_cdf_g <- pgamma(mass_seq, shape=mass_g$estimate[1],
rate=mass_g$estimate[2])
mass_cdf <- ecdf(Fl_data$mass)(mass_seq)
mass_cdf_df <- as.data.frame(mass_seq) %>%
mutate(cdf_ln = mass_cdf_ln,
cdf_n = mass_cdf_n,
cdf_g = mass_cdf_g,
cdf = mass_cdf)
## Plot Mass PDF
ggplot(mass_pdf_df)+
geom_point(aes(x=mass_seq, y=pdf_n, color="Normal"))+
geom_point(aes(x=mass_seq, y=pdf_ln, color="Lognormal"))+
geom_point(aes(x=mass_seq, y=pdf_g, color="Gamma"))+
geom_point(aes(x=mass_pdf$x, y=mass_pdf$y, color="Data"))+
theme_bw()+
scale_color_manual(values=four_colors)+
labs(title="PDF for mass data and possible distributions")+
ylab("PDF")+
xlab("mass")+
theme(
text = element_text(size=24, family="mont"),
legend.position="inside",
legend.position.inside = c(.85,.6),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
## Plot Mass CDF
ggplot(mass_cdf_df)+
geom_point(aes(x=mass_seq, y=cdf_n, color="Normal"))+
geom_point(aes(x=mass_seq, y=cdf_ln, color="Lognormal"))+
geom_point(aes(x=mass_seq, y=cdf_g, color="Gamma"))+
geom_point(aes(x=mass_seq, y=cdf, color="Data"))+
theme_bw()+
scale_color_manual(values=four_colors)+
labs(title="CDF for mass data and possible distributions")+
ylab("CDF")+
xlab("mass")+
theme(
text = element_text(size=24, family="mont"),
legend.position="inside",
legend.position.inside = c(.85,.6),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
### test mass for gamma and plnorm distributions using ks test
ks.test(Fl_data$mass, "pgamma",shape=mass_g$estimate[1],
rate=mass_g$estimate[2])## Warning in ks.test.default(Fl_data$mass, "pgamma", shape = mass_g$estimate[1],
## : ties should not be present for the Kolmogorov-Smirnov test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: Fl_data$mass
## D = 0.05821, p-value = 0.001924
## alternative hypothesis: two-sidedks.test(Fl_data$mass, "plnorm",meanlog=mass_ln$estimate[1],
sdlog=mass_ln$estimate[2])## Warning in ks.test.default(Fl_data$mass, "plnorm", meanlog =
## mass_ln$estimate[1], : ties should not be present for the Kolmogorov-Smirnov
## test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: Fl_data$mass
## D = 0.039922, p-value = 0.07622
## alternative hypothesis: two-sided## sugar_avg distribution
sugar_avg_n <- fitdistr(Fl_data_no_BER$sugar_avg, "normal")
sugar_avg_ln <- fitdistr(Fl_data_no_BER$sugar_avg, "lognormal")
sugar_avg_g <- fitdistr(Fl_data_no_BER$sugar_avg, "gamma")## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced## Warning in densfun(x, parm[1], parm[2], ...): NaNs produced
## Warning in densfun(x, parm[1], parm[2], ...): NaNs producedsugar_avg_seq <- seq(min(Fl_data_no_BER$sugar_avg), max(Fl_data_no_BER$sugar_avg), by=0.01)
## PDFs
sugar_avg_pdf_n <- dnorm(sugar_avg_seq, mean=sugar_avg_n$estimate[1],
sd=sugar_avg_n$estimate[2])
sugar_avg_pdf_ln <- dlnorm(sugar_avg_seq, meanlog=sugar_avg_ln$estimate[1],
sdlog=sugar_avg_ln$estimate[2])
sugar_avg_pdf_g <- dgamma(sugar_avg_seq, shape=sugar_avg_g$estimate[1],
rate=sugar_avg_g$estimate[2])
sugar_avg_pdf <- density(Fl_data_no_BER$sugar_avg, n=length(sugar_avg_pdf_ln))
sugar_avg_pdf_df <- as.data.frame(sugar_avg_seq) %>%
mutate(pdf_ln = sugar_avg_pdf_ln,
pdf_n = sugar_avg_pdf_n,
pdf_g = sugar_avg_pdf_g)
## CDFs
sugar_avg_cdf_n <- pnorm(sugar_avg_seq, mean=sugar_avg_n$estimate[1],
sd=sugar_avg_n$estimate[2])
sugar_avg_cdf_ln <- plnorm(sugar_avg_seq, meanlog=sugar_avg_ln$estimate[1],
sdlog=sugar_avg_ln$estimate[2])
sugar_avg_cdf_g <- pgamma(sugar_avg_seq, shape=sugar_avg_g$estimate[1],
rate=sugar_avg_g$estimate[2])
sugar_avg_cdf <- ecdf(Fl_data$sugar_avg)(sugar_avg_seq)
sugar_avg_cdf_df <- as.data.frame(sugar_avg_seq) %>%
mutate(cdf_ln = sugar_avg_cdf_ln,
cdf_n = sugar_avg_cdf_n,
cdf_g = sugar_avg_cdf_g,
cdf = sugar_avg_cdf)
## Plot sugar PDF
ggplot(sugar_avg_pdf_df)+
geom_point(aes(x=sugar_avg_seq, y=pdf_n, color="Normal"))+
geom_point(aes(x=sugar_avg_seq, y=pdf_ln, color="Lognormal"))+
geom_point(aes(x=sugar_avg_seq, y=pdf_g, color="Gamma"))+
geom_point(aes(x=sugar_avg_pdf$x, y=sugar_avg_pdf$y, color="Data"))+
theme_bw()+
scale_color_manual(values=four_colors)+
labs(title="PDF for sugar_avg data and possible distributions")+
ylab("PDF")+
xlab("sugar_avg")+
theme(
text = element_text(size=24, family="mont"),
legend.position="inside",
legend.position.inside = c(.85,.6),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
## Plot sugar CDF
ggplot(sugar_avg_cdf_df)+
geom_point(aes(x=sugar_avg_seq, y=cdf_n, color="Normal"))+
geom_point(aes(x=sugar_avg_seq, y=cdf_ln, color="Lognormal"))+
geom_point(aes(x=sugar_avg_seq, y=cdf_g, color="Gamma"))+
geom_point(aes(x=sugar_avg_seq, y=cdf, color="Data"))+
theme_bw()+
scale_color_manual(values=four_colors)+
labs(title="CDF for sugar_avg data and possible distributions")+
ylab("CDF")+
xlab("sugar_avg")+
theme(
text = element_text(size=24, family="mont"),
legend.position="inside",
legend.position.inside = c(.85,.6),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold", lineheight = .5)
)
### test sugar_avg for gamma and plnorm distributions using ks test
ks.test(Fl_data_no_BER$sugar_avg, "pgamma",shape=sugar_avg_g$estimate[1],
rate=sugar_avg_g$estimate[2])## Warning in ks.test.default(Fl_data_no_BER$sugar_avg, "pgamma", shape =
## sugar_avg_g$estimate[1], : ties should not be present for the
## Kolmogorov-Smirnov test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: Fl_data_no_BER$sugar_avg
## D = 0.068149, p-value = 0.005857
## alternative hypothesis: two-sidedks.test(Fl_data_no_BER$sugar_avg, "plnorm",meanlog=sugar_avg_ln$estimate[1],
sdlog=sugar_avg_ln$estimate[2])## Warning in ks.test.default(Fl_data_no_BER$sugar_avg, "plnorm", meanlog =
## sugar_avg_ln$estimate[1], : ties should not be present for the
## Kolmogorov-Smirnov test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: Fl_data_no_BER$sugar_avg
## D = 0.083045, p-value = 0.0003461
## alternative hypothesis: two-sidedks.test(Fl_data_no_BER$sugar_avg, "pnorm",mean=sugar_avg_n$estimate[1],
sd=sugar_avg_n$estimate[2])## Warning in ks.test.default(Fl_data_no_BER$sugar_avg, "pnorm", mean =
## sugar_avg_n$estimate[1], : ties should not be present for the
## Kolmogorov-Smirnov test##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: Fl_data_no_BER$sugar_avg
## D = 0.062847, p-value = 0.01401
## alternative hypothesis: two-sided## Test for homoscedasticity
leveneTest(Fl_data_no_BER$sugar_avg~Fl_data_no_BER$Treatment)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 2.211 0.08568 .
## 624
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1leveneTest(Fl_data_no_BER$mass~Fl_data_no_BER$Treatment)## Levene's Test for Homogeneity of Variance (center = median)
## Df F value Pr(>F)
## group 3 1.7838 0.149
## 624# Plots
## Density plot for mass by treatment across no BER tomatoes (exploratory)
ggplot(data=Fl_data_no_BER, aes(x=mass, y = Treatment, color=Treatment, fill=Treatment))+
geom_density_ridges()+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
#xlim(-0.25, 1.5)+
ylab("Density")+
xlab("Mass(g)")+
theme_minimal()+
labs(
title=str_wrap("Mass Density for Inoculation Treatments in Tomato", 40)
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight = 0.5),
axis.title.y = element_text(size=0),
axis.text = element_text(size=24, family="mont"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Picking joint bandwidth of 23.2
## Density plot for mass by plant across no BER tomatoes (exploratory)
ggplot(data=Fl_data_no_BER, aes(x=mass, y = plant_fac, color=Treatment, fill=Treatment))+
geom_density_ridges()+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
xlim(0, 400)+
ylab("Plant")+
xlab("Mass (g)")+
theme_minimal()+
labs(
title=str_wrap("Mass Density for Individual Plants Across Inoculation Treatments in Tomato", 50)
)+
theme(
legend.position="right",
text = element_text(size=24, family="mont", lineheight = 0.5),
axis.text.x = element_text(size=24, family="mont"),
axis.text.y = element_text(size=12, family="mont"),
axis.title = element_text(size=28, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Picking joint bandwidth of 30.9## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_density_ridges()`).
## Density plots for sugar by treatment across no BER tomatoes (exploratory)
ggplot(data=Fl_data_no_BER, aes(x=sugar_avg, y = Treatment, color=Treatment, fill=Treatment))+
geom_density_ridges()+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
#xlim(-0.25, 1.5)+
ylab("Density")+
xlab("Sugar Concentration (%)")+
theme_minimal()+
labs(
title=str_wrap("Sugar Concentration Density for Inoculation Treatments in Tomato", 40)
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", lineheight = 0.5),
axis.title.y = element_text(size=0),
axis.text = element_text(size=24, family="mont"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Picking joint bandwidth of 0.388
## Density plots for sugar by plant across no BER tomatoes (exploratory)
ggplot(data=Fl_data_no_BER, aes(x=sugar_avg, y = plant_fac, color=Treatment, fill=Treatment))+
geom_density_ridges()+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
xlim(3, 9)+
ylab("Plant")+
xlab("Sugar Concentration (%)")+
theme_minimal()+
labs(
title=str_wrap("Sugar Concentration Density for Individual Plants Across Inoculation Treatments in Tomato", 50)
)+
theme(
legend.position="right",
text = element_text(size=24, family="mont", lineheight = 0.5),
axis.text.y = element_text(size=12, family="mont"),
axis.title = element_text(size=28, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Picking joint bandwidth of 0.447## Warning: Removed 51 rows containing non-finite outside the scale range
## (`stat_density_ridges()`).
## scatter plot of sugar by mass grouped and facet wrapped by treatment
ggplot(data=Fl_data_no_BER, aes(x=mass, y=sugar_avg,
fill=Treatment, shape=Treatment,
color=Treatment))+
geom_point()+
scale_color_manual(values=four_colors)+
scale_fill_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
ylab("Sugar Concentration (%)")+
xlab("Tomato Mass (g)")+
theme_minimal()+
facet_wrap(~Treatment)+
labs(
title="Tomato Sugar Concentration by Mass"
)+
theme(
legend.position="right",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## scatter plot of ripeness by 'days from harvest to analysis' grouped by treatment (exploratory, to see how delayed analysis may affect sugar content [hopefully it doesn't])
ggplot(data=Fl_data_no_BER, aes(x=d_diff, y=ripeness,
fill=Treatment, shape=Treatment,
color=Treatment))+
geom_jitter()+
scale_color_manual(values=four_colors)+
scale_fill_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
ylab("Ripeness")+
xlab("Days from Harvest to Analysis")+
theme_minimal()+
#facet_wrap(~Treatment)+
labs(
title="Tomato Ripeness by Days from Harvest to Analysis"
)+
theme(
legend.position="right",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Inf## Warning in min(d[d > tolerance]): no non-missing arguments to min; returning
## Inf## Warning in min(x): no non-missing arguments to min; returning Inf## Warning in max(x): no non-missing arguments to max; returning -Inf## Warning in stats::runif(length(x), -amount, amount): NAs produced## Warning: Removed 628 rows containing missing values or values outside the scale range
## (`geom_point()`).
## scatter plot of sugar by days from harvest to analysis grouped by treatment (exploratory, to see how delayed analysis may affect sugar content [hopefully it doesn't])
ggplot(data=Fl_data_no_BER, aes(x=d_diff, y=sugar_avg,
fill=Treatment, shape=Treatment,
color=Treatment))+
geom_jitter()+
scale_fill_manual(values=four_colors)+
scale_shape_manual(values=four_shapes)+
scale_color_manual(values=four_colors)+
ylab("Sugar Concentration (%)")+
xlab("Days from Harvest to Analysis")+
theme_minimal()+
#facet_wrap(~Treatment)+
labs(
title="Tomato Sugar Concentration by Days from Harvest to Analysis"
)+
theme(
legend.position="right",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
3.2.2 Fruit Occurrence Analysis
Looking at the occurrence of BER, fungus, and cracking by various groups.
Because of how few counts we had for certain factors (<50 observations for each group in BER, cracking, fungus), a chi-square was not suitable, so Fisher’s exact test was used.
## Treatment grouped fisher test and row-wise fisher test for BER
sum_BER_tab <- as.table(cbind(Fl_data_summary$Treatment,
Fl_data_summary$BER_sum))
fisher.test(Fl_data_summary$BER_sum, y=Fl_data_summary$Treatment)##
## Fisher's Exact Test for Count Data
##
## data: Fl_data_summary$BER_sum and Fl_data_summary$Treatment
## p-value = 1
## alternative hypothesis: two.sidedrow_wise_fisher_test(xtab = sum_BER_tab)## # A tibble: 4 × 5
## group n p p.adj p.adj.signif
## * <chr> <int> <dbl> <dbl> <chr>
## 1 A 120 0.688 1 ns
## 2 B 120 0.721 1 ns
## 3 C 120 0.725 1 ns
## 4 D 120 0.221 0.884 ns## Treatment grouped fishers exact test on fruit
sum_fruit_tab <- as.table(cbind(Fl_data_summary$Treatment,
Fl_data_summary$fruit_sum))
fisher.test(Fl_data_summary$fruit_sum, y=Fl_data_summary$Treatment)##
## Fisher's Exact Test for Count Data
##
## data: Fl_data_summary$fruit_sum and Fl_data_summary$Treatment
## p-value = 1
## alternative hypothesis: two.sidedrow_wise_fisher_test(xtab = sum_fruit_tab)## # A tibble: 4 × 5
## group n p p.adj p.adj.signif
## * <chr> <dbl> <dbl> <dbl> <chr>
## 1 A 1035 0.697 1 ns
## 2 B 1035 1 1 ns
## 3 C 1035 1 1 ns
## 4 D 1035 0.122 0.488 ns## Treatment grouped fishers exact test on fungus
sum_fungus_tab <- as.table(cbind(Fl_data_summary$Treatment,
Fl_data_summary$fungus_sum))
fisher.test(Fl_data_summary$fungus_sum, y=Fl_data_summary$Treatment)##
## Fisher's Exact Test for Count Data
##
## data: Fl_data_summary$fungus_sum and Fl_data_summary$Treatment
## p-value = 1
## alternative hypothesis: two.sidedrow_wise_fisher_test(xtab = sum_fungus_tab)## # A tibble: 4 × 5
## group n p p.adj p.adj.signif
## * <chr> <int> <dbl> <dbl> <chr>
## 1 A 187 0.287 0.861 ns
## 2 B 187 1 1 ns
## 3 C 187 1 1 ns
## 4 D 187 0.102 0.408 ns## Treatment grouped fishers exact test on cracking
sum_cracking_tab <- as.table(cbind(Fl_data_summary$Treatment,
Fl_data_summary$cracking_sum))
fisher.test(Fl_data_summary$cracking_sum, y=Fl_data_summary$Treatment)##
## Fisher's Exact Test for Count Data
##
## data: Fl_data_summary$cracking_sum and Fl_data_summary$Treatment
## p-value = 1
## alternative hypothesis: two.sidedrow_wise_fisher_test(xtab = sum_cracking_tab)## # A tibble: 4 × 5
## group n p p.adj p.adj.signif
## * <chr> <int> <dbl> <dbl> <chr>
## 1 A 273 0.462 1 ns
## 2 B 273 1 1 ns
## 3 C 273 0.749 1 ns
## 4 D 273 0.225 0.9 ns## Mass bin grouped fishers exact test on BER
mb_BER_tab <- as.table(cbind(Fl_data_mb$Treatment, Fl_data_mb$BER_sum))
fisher.test(Fl_data_mb$BER_sum, y=Fl_data_mb$Treatment)##
## Fisher's Exact Test for Count Data
##
## data: Fl_data_mb$BER_sum and Fl_data_mb$Treatment
## p-value = 0.9962
## alternative hypothesis: two.sidedrow_wise_fisher_test(xtab = mb_BER_tab)## # A tibble: 37 × 5
## group n p p.adj p.adj.signif
## * <chr> <int> <dbl> <dbl> <chr>
## 1 A 202 0.223 1 ns
## 2 B 202 0.129 1 ns
## 3 C 202 0.0071 0.256 ns
## 4 D 202 1 1 ns
## 5 E 202 1 1 ns
## 6 F 202 0.455 1 ns
## 7 G 202 0.455 1 ns
## 8 H 202 1 1 ns
## 9 I 202 0.455 1 ns
## 10 J 202 0.0000135 0.000499 ***
## # ℹ 27 more rows## Mass bin grouped fishers exact test on fruit
mb_fruit_tab <- as.table(cbind(Fl_data_mb$Treatment, Fl_data_mb$fruit_sum))
fisher.test(Fl_data_mb$fruit_sum, y=Fl_data_mb$Treatment)##
## Fisher's Exact Test for Count Data
##
## data: Fl_data_mb$fruit_sum and Fl_data_mb$Treatment
## p-value = 0.5413
## alternative hypothesis: two.sidedrow_wise_fisher_test(xtab = mb_fruit_tab)## # A tibble: 37 × 5
## group n p p.adj p.adj.signif
## * <chr> <dbl> <dbl> <dbl> <chr>
## 1 A 1117 0.505 1 ns
## 2 B 1117 0.118 1 ns
## 3 C 1117 0.0811 1 ns
## 4 D 1117 0.715 1 ns
## 5 E 1117 1 1 ns
## 6 F 1117 1 1 ns
## 7 G 1117 0.404 1 ns
## 8 H 1117 0.291 1 ns
## 9 I 1117 0.158 1 ns
## 10 J 1117 0.0025 0.0825 ns
## # ℹ 27 more rows## Mass bin grouped fishers exact test on fungus
mb_fungus_tab <- as.table(cbind(Fl_data_mb$Treatment, Fl_data_mb$fungus_sum))
fisher.test(Fl_data_mb$fungus_sum, y=Fl_data_mb$Treatment)##
## Fisher's Exact Test for Count Data
##
## data: Fl_data_mb$fungus_sum and Fl_data_mb$Treatment
## p-value = 0.8642
## alternative hypothesis: two.sidedrow_wise_fisher_test(xtab = mb_fungus_tab)## # A tibble: 37 × 5
## group n p p.adj p.adj.signif
## * <chr> <int> <dbl> <dbl> <chr>
## 1 A 269 0.172 1 ns
## 2 B 269 0.0396 1 ns
## 3 C 269 0.00164 0.0607 ns
## 4 D 269 0.667 1 ns
## 5 E 269 0.172 1 ns
## 6 F 269 1 1 ns
## 7 G 269 0.342 1 ns
## 8 H 269 0.342 1 ns
## 9 I 269 0.342 1 ns
## 10 J 269 0.0242 0.774 ns
## # ℹ 27 more rows## Mass bin grouped fishers exact test on cracking
mb_cracking_tab <- as.table(cbind(Fl_data_mb$Treatment,
Fl_data_mb$cracking_sum))
fisher.test(Fl_data_mb$cracking_sum, y=Fl_data_mb$Treatment)##
## Fisher's Exact Test for Count Data
##
## data: Fl_data_mb$cracking_sum and Fl_data_mb$Treatment
## p-value = 0.8969
## alternative hypothesis: two.sidedrow_wise_fisher_test(xtab = mb_cracking_tab)## # A tibble: 37 × 5
## group n p p.adj p.adj.signif
## * <chr> <int> <dbl> <dbl> <chr>
## 1 A 355 1 1 ns
## 2 B 355 0.0512 1 ns
## 3 C 355 0.0201 0.663 ns
## 4 D 355 0.686 1 ns
## 5 E 355 0.456 1 ns
## 6 F 355 1 1 ns
## 7 G 355 1 1 ns
## 8 H 355 0.259 1 ns
## 9 I 355 0.259 1 ns
## 10 J 355 0.258 1 ns
## # ℹ 27 more rows## probability of BER for each mass bin
ggplot(data=Fl_data_mb, aes(x=mass_bin, y=pBER, fill=mass_bin))+
geom_col()+
ylim(0, 50)+
scale_fill_manual(values=ten_col)+
ylab("BER Probability (%)")+
xlab("Mass Bin (g)")+
theme_minimal()+
labs(
title="Tomato BER Probability by Mass Bin"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_col()`).
## Plot probability of blossom end rot by treatment as a column graph
ggplot(data=Fl_data_summary, aes(x=Treatment, y=pBER*100, fill=Treatment))+
geom_col()+
scale_fill_manual(values=four_colors)+
ylab("Percent of Fruit w/ BER")+
xlab("Treatment")+
theme_minimal()+
labs(
title="Tomato Fruit Blossom End Rot Occurrence by Treatment"
)+
theme(
legend.position="none",
text = element_text(size=30, family="mont", face="bold"),
axis.title = element_text(size=30, family = "mont", face= "bold"),
title = element_text(size=34, family="open", face="bold")
)
## Plot probability of blossom end rot by treatment as a column graph
ggplot(data=Fl_data_summary, aes(x=Treatment, y=fruit_sum, fill=Treatment))+
geom_col()+
scale_fill_manual(values=four_colors)+
ylab("Fruit Count")+
xlab("Treatment")+
theme_minimal()+
labs(
title="Total Fruit Count by Inoculation Treatment"
)+
theme(
legend.position="none",
text = element_text(size=30, family="mont", face="bold"),
axis.title = element_text(size=30, family = "mont", face= "bold"),
title = element_text(size=34, family="open", face="bold")
)
## Plot blossom end rot by treatment as a stacked bar chart
ggplot(data=Fl_data, aes(x=Treatment, y=Fl_data[order(Fl_data$BER),]$fruit, fill=BER_fac))+
geom_bar(position="stack", stat="identity")+
scale_fill_manual(values=two_colors)+
ylab("Fruit Count")+
xlab("Treatment")+
theme_minimal()+
guides(fill=guide_legend(title=str_wrap("Blossom end rot", 8)))+
labs(
title="Tomato Fruit and BER Count by Treatment"
)+
theme(
legend.position="right",
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## Plot blossom end rot by mass bin as a stacked bar chart
ggplot(data=Fl_data, aes(x=mass_bin, y=Fl_data[order(Fl_data$BER),]$fruit, fill=BER_fac))+
geom_bar(position="stack", stat="identity")+
scale_fill_manual(labels=c("No", "Yes"), values=two_colors)+
guides(fill=guide_legend(title=str_wrap("Blossom end rot", 8)))+
ylab("Fruit Count")+
xlab("Mass Bin")+
theme_minimal()+
labs(
title="Tomato Fruit and BER Count by Mass Bin"
)+
theme(
legend.position="inside",
legend.position.inside=c(.8,.7),
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
axis.text.x = element_text(size=20, family = "mont", angle=45,
lineheight=1, hjust=1, vjust=1.3),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## Plot blossom end rot by plant as a stacked bar chart
ggplot(data=Fl_data, aes(x=plant, y=Fl_data[order(Fl_data$BER),]$fruit, fill=BER_fac))+
geom_bar(position="stack", stat="identity")+
scale_fill_manual(labels=c("No", "Yes"), values=two_colors)+
guides(fill=guide_legend(title=str_wrap("Blossom end rot", 8)))+
ylab("Fruit Count")+
xlab("Plant")+
facet_wrap(~Treatment)+
theme_minimal()+
labs(
title="Tomato Fruit and BER Count by Plant"
)+
theme(
legend.position="right",
legend.background = element_rect(color=four_colors[3], fill=NA,
linewidth=.5, linetype = 2),
text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
axis.text = element_text(size=16, family = "mont",
lineheight=1),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## fruit mass sum by treatment
ggplot(data=Fl_data_summary, aes(x=Treatment, y=mass_sum/1000,
fill=Treatment))+
geom_col()+
scale_fill_manual(values=four_colors)+
ylab("Total Fruit Mass (kg)")+
xlab("Treatment")+
theme_minimal()+
labs(
title="Total Tomato Fruit Harvest across Inoculation Treatments"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## fruit mass sum by plant
ggplot(data=Fl_sum_byplant, aes(x=plant, y=mass_sum/1000,
fill=Treatment))+
geom_col()+
facet_wrap(~Treatment)+
scale_fill_manual(values=four_colors)+
scale_x_discrete(name="Plant", breaks=c("3","6","9"))+
ylab("Total Fruit Mass (kg)")+
xlab("Plant")+
theme_minimal()+
labs(
title="Total Tomato Fruit Harvest across Individual Plants"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## plot of fruit mass by treatment
ggplot(data=Fl_sum_byplant, aes(x=Treatment, y=mass_sum/1000, color=Treatment, fill=Treatment))+
geom_violin(alpha=0.5)+
geom_boxplot(width=0.2, alpha=0.5)+
geom_jitter(width=0.1)+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Mass Sum (kg)")+
xlab("Treatment")+
theme_minimal()+
labs(
title="Total Fruit Mass per Plant across Inoculation Treatments"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## plot of fruit count by treatment
ggplot(data=Fl_sum_byplant, aes(x=Treatment, y=fruit_sum, color=Treatment, fill=Treatment))+
geom_violin(alpha=0.5)+
geom_boxplot(width=0.2, alpha=0.5)+
geom_jitter(width=0.1)+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Fruit Count")+
xlab("Treatment")+
theme_minimal()+
labs(
title="Tomato Plant Harvested Fruit Count Across Inoculation Treatments"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
Fl_avg_sums <- Fl_sum_byplant %>%
group_by(Treatment) %>%
summarize_at(vars(mass_sum), list(mean=mean, sd=sd))%>%
mutate(mean=mean/1000,
sd=sd/1000)
## mass per treatment columns with error bars
ggplot(data=Fl_avg_sums, aes(x=Treatment, y=mean, fill=Treatment))+
geom_col()+
geom_errorbar( aes(ymin=mean-sd, ymax=mean+sd), width=0.4)+
scale_fill_manual(values=four_colors)+
ylab("Total Mass (kg)")+
xlab("Treatment")+
theme_minimal()+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
3.2.3 Mass Analysis
#Kruskal-Wallis on mass as a function of Treatment
kruskal.test(mass ~ Treatment, data=Fl_data)##
## Kruskal-Wallis rank sum test
##
## data: mass by Treatment
## Kruskal-Wallis chi-squared = 57.667, df = 3, p-value = 1.851e-12kruskal.test(mass_sum ~ Treatment, data=Fl_sum_byplant)##
## Kruskal-Wallis rank sum test
##
## data: mass_sum by Treatment
## Kruskal-Wallis chi-squared = 9.5145, df = 3, p-value = 0.02318# Pairwise wilcox test on mass as a function of treatment
pairwise.wilcox.test(Fl_data$mass, Fl_data$Treatment)##
## Pairwise comparisons using Wilcoxon rank sum test with continuity correction
##
## data: Fl_data$mass and Fl_data$Treatment
##
## Control Germination Transplantation
## Germination 3.8e-10 - -
## Transplantation 3.8e-10 0.94819 -
## Germ+Trans 0.00019 0.00807 0.00807
##
## P value adjustment method: holmpairwise.wilcox.test(Fl_sum_byplant$mass_sum, Fl_sum_byplant$Treatment)##
## Pairwise comparisons using Wilcoxon rank sum exact test
##
## data: Fl_sum_byplant$mass_sum and Fl_sum_byplant$Treatment
##
## Control Germination Transplantation
## Germination 1.000 - -
## Transplantation 0.239 0.086 -
## Germ+Trans 1.000 1.000 0.014
##
## P value adjustment method: holmmT_glm_s <- summary(mT_glm <- glm(mass~Treatment, data=Fl_data, family=Gamma(link="log")))
mT_glm_pseudoR2 <- (mT_glm_s$null.deviance - mT_glm_s$deviance)/mT_glm_s$null.deviance
summary(glm(mass_sum~Treatment, data=Fl_sum_byplant, family=Gamma(link="log")))##
## Call:
## glm(formula = mass_sum ~ Treatment, family = Gamma(link = "log"),
## data = Fl_sum_byplant)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.52975 0.10213 73.727 <2e-16 ***
## TreatmentGermination -0.03063 0.14443 -0.212 0.8330
## TreatmentTransplantation 0.29324 0.14443 2.030 0.0484 *
## TreatmentGerm+Trans -0.15248 0.14443 -1.056 0.2969
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for Gamma family taken to be 0.1251669)
##
## Null deviance: 7.0470 on 47 degrees of freedom
## Residual deviance: 5.7134 on 44 degrees of freedom
## AIC: 764.79
##
## Number of Fisher Scoring iterations: 4## plot of mass by treatment with box plots, violins, and jittered points
ggplot(data=Fl_data_no_BER, aes(x=Treatment, y=mass, color=Treatment, fill=Treatment))+
geom_violin(alpha=0.5)+
geom_boxplot(width=0.2, alpha=0.5)+
geom_jitter(width=0.1)+
stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
stat_compare_means(label.x=3.25, label.y=425, size=8, family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Tomato Mass (g)")+
xlab("Treatment")+
theme_minimal()+
#facet_wrap(~Treatment)+
labs(
title="Tomato Mass by Microbial Inoculation Timing"
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold"),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
## plot of mass mean per plant by treatment with box plots, violins, and jittered points
ggplot(data=Fl_means_byplant, aes(x=Treatment, y=mass_mean, color=Treatment, fill=Treatment))+
geom_violin(alpha=0.5)+
geom_boxplot(width=0.2, alpha=0.5)+
geom_jitter(width=0.1)+
stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
stat_compare_means(label.x=3.25, size=8, family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Tomato Mass (g)")+
xlab("Treatment")+
theme_minimal()+
#facet_wrap(~Treatment)+
labs(
title=str_wrap("Mean Tomato Mass by Plant across Microbial Inoculation Timing", 50)
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
3.2.4 Sugar Analysis
# Linear mixed model of sugar as a function of treatment with plant as a random effect
#summary(model_sugar_no_BER <- (lmer(
# sugar_avg ~ Treatment + (1| plant_fac),
# data = Fl_data_no_BER)))
#summary(lm_sugar_no_BER <- (lm(
# sugar_avg ~ Treatment,
# data = Fl_data_no_BER)))
#TUKEY (more or less) on the LMER
#difflsmeans(model_sugar_no_BER)
## plot sugar by treatment as a boxplot with violins and jitter
ggplot(data=Fl_data_no_BER, aes(x=Treatment, y=sugar_avg, color=Treatment, fill=Treatment))+
geom_violin(alpha=0.5)+
geom_boxplot(width=0.2, alpha=0.5)+
geom_jitter(width=0.1)+
stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
stat_compare_means(label.x=3.5, label.y=11, size=8, family="mont")+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Sugar Concentration (%)")+
xlab("Treatment")+
theme_minimal()+
#facet_wrap(~Treatment)+
labs(
title=str_wrap("Tomato Sugar Concentration by Microbial Inoculation Timing", 40)
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)
sugar_kw <- kruskal.test(sugar_avg_mean~Treatment, Fl_means_byplant)
## plot sugar means per plant by treatment as a boxplot with violins and jitter
ggplot(data=Fl_means_byplant, aes(x=Treatment, y=sugar_avg_mean, color=Treatment, fill=Treatment))+
geom_violin(alpha=0.5)+
geom_boxplot(width=0.2, alpha=0.5)+
geom_jitter(width=0.1)+
stat_compare_means(comparisons = list(c("Control","Germination"),c("Control", "Transplantation"), c("Control", "Germ+Trans")), size=8, family="mont")+
annotate("text", x=3.5, y=8.5, family="mont", label= paste("Kruskal-Wallis, p=", round(sugar_kw$p.value, 4)), size=8)+
scale_fill_manual(values=four_colors)+
scale_color_manual(values=four_colors)+
ylab("Sugar Concentration (%)")+
xlab("Treatment")+
theme_minimal()+
#facet_wrap(~Treatment)+
labs(
title=str_wrap("Tomato Plant Sugar Concentration Mean by Microbial Inoculation Timing", 40)
)+
theme(
legend.position="none",
text = element_text(size=24, family="mont", face="bold", lineheight=0.5),
axis.title = element_text(size=24, family = "mont", face= "bold"),
title = element_text(size=30, family="open", face="bold")
)